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Abstract 

Phenotypic fluctuations and plasticity can generally affect the course 
of evolution, a process known as the Baldwin effect. Several studies 
have recast this effect and claimed that phenotypic plasticity acceler- 
ates evolutionary rate (the Baldwin expediting effect); however, the 
validity of this claim is still controversial. In this study, we investi- 
gate the evolutionary population dynamics of a quantitative genetic 
model under a multi-peaked fitness landscape, in order to evaluate 
the validity of the effect. We provide analytical expressions for the 
evolutionary rate and average population fitness. Our results indicate 
that under a multi-peaked fitness landscape, phenotypic fluctuation 
always accelerates evolutionary rate, but it decreases the average fit- 
ness. As an extreme case of the trade-off between the rate of evolution 
and average fitness, phenotypic fiuctuation is shown to accelerate the 
error catastrophe, in which a population fails to sustain a high-fitness 
peak. In the context of our findings, we discuss the role of phenotypic 
plasticity in adaptive evolution. 
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1 Introduction 



Phenotypic plasticity is a common process, whereby phenotypic changes oc- 
cur during an organism's hfe cycle. In some cases, phenotypic plasticity can 
emerge as a response to environmental variation. For example, the desert 
locust transitions from a solitary to a gregarious phase, depending on the 
density of population [ll[2]. In addition, a phenotype can also fluctuate ran- 
domly, for example, as a result of stochastic nature in gene expression, as 
exemplified by recent studies in microorganisms [31 HI [5], [6] and animal devel- 
opment [ZllHllQ]. Such processes are also interpreted as phenotypic plasticity, 
and are ubiquitous in nature. Given this, the role of random phenotypic 
plasticity in adaptive evolution warrants further consideration. 

Do phenotypic changes in individuals have any impacts over generations, 
and thus, over evolution? At first glance, phenotypic plasticity does not seem 
to affect evolutionary processes because only genotypes, rather than pheno- 
types, are heritable. However, non-heritable phenotypes acquired during the 
life cycle can be fixed into genotypes through natural selection acting on 
phenotypes, and thus, can affect evolution. This process is referred to as the 
Baldwin effect, which was first introduced by Baldwin [10], and later recast 
by Simpson [11]. Several studies [12], [131 [HI [15], O [161 [17] have revealed that 
this effect is advantageous under a fluctuating environment, as phenotypic 
plasticity can allow organisms to survive sudden environmental changes, and 
therefore, facilitate adaptation to the new environment. However, whether 
phenotypic plasticity is advantageous throughout evolution, even under a 
fixed environment, is still elusive. 

Can phenotypic plasticity accelerate evolutionary rate? This question 
stems from a recent incarnation of the Baldwin effect known as the " Baldwin 
expediting effect" [18]. The seminal work by Hinton and Nowlan [19], indeed, 
showed that phenotypic plasticity accelerates evolutionary rate even under 
a fixed environment, although, the validity of these findings is controversial. 
For example, several subsequent studies have shown that phenotypic plastic- 
ity can decelerate evolutionary rate [Sni [El [21] , whereas others have claimed 
that it accelerates evolutionary rate [19], [22] [23] [211 ES, [26]. However, it is 
important to note that the studies in which it was concluded that phenotypic 
plasticity leads to a deceleration of evolutionary rate dealt with uni-modal 
fitness landscapes. Recently, some studies pointed out that, under a multi- 
peaked fitness, phenotypic plasticity can smoothen fitness valleys, and thus, 
the escape from a local maximum can be enhanced, leading to the accelera- 
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tion of evolution [271 ESI [2H1 ES]- In fact, the acceleration of evolution under 
a multi-peaked fitness landscape has been confirmed numerically [221 [26] . 

The first analytical result addressing evolution under a multi-peaked fit- 
ness landscape in a fixed environment was provided by Borenstein et al. [25] . 
in which they showed that phenotypic plasticity accelerates evolution. In- 
stead of a mutation-selection process, they considered a single random walker 
on a multi-peaked fitness landscape, to mimic evolutionary population dy- 
namics in genotypic space; therefore, their study did not distinguish between 
the behaviors of a population and an individual. Thus in their random walker 
model the following three points were left unaddressed; first, average fitness 
over population has not been obtained analytically; second, behaviors of the 
distribution of population cannot be discussed in their model; third, the rela- 
tion between the acceleration of evolution and average fitness over population 
is unclear. These points should be clarified to fully comprehend the role of 
phenotypic fluctuation (plasticity) in evolutionary population dynamics. 

In this study, we investigate evolutionary population dynamics, rather 
than the single random walker model, under multi-peaked landscapes in a 
fixed environment and study the role of phenotypic fluctuation on the evo- 
lutionary rate. We provide analytical results for the relationship between 
phenotypic fluctuation, the evolutionary rate, and the average population 
fitness by considering the distribution of genotypes in a population. First, 
we assess whether a population located at a fitness peak can stay around 
the peak over generations to maintain high fitness. We show that for both 
large phenotypic fluctuations and high mutation rates, populations fail to 
keep descendants near the peak over multiple generations, and thus are un- 
able to maintain a high average fitness. Second, we evaluate the rate at 
which a population moves between fitness peaks, and derive an analytical 
expression representing the dependency of evolutionary rate on phenotypic 
fluctuation. The obtained result indicates that phenotypic fluctuation always 
enhances diffusion in genotypic space, and thus accelerates evolution under 
multi-peaked fitness landscapes. 

One consequence of our results is that, for a large phenotypic fluctua- 
tion, the failure of a population to concentrate near a fintness peak results 
in a decrease in average fitness. This can be regarded as "error catastro- 
phe." The concept of error catastrophe refers to the collapse of a population 
that maintains high fitness because of a high mutation rate [301 EI]. Our 
results demonstrate that the error catastrophe occurs not only because of 
a high mutation rate, but also as a result of large phenotypic fluctuations. 
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Furthermore, our results illustrate an apparent relationship between error 
catastrophe and the Baldwin effect, in that error catastrophe is explained as 
an extreme case of a trade-off between an acceleration of evolutionary rate 
and a decrease in fitness. These findings reveal novel aspects of the role of 
phenotypic plasticity in evolution. 

The organization of this paper is as follows. In Sec. 12.11 we explain how 
phenotypic fluctuation modifies fitness landscapes. In Sec. 12.21 we introduce 
our model of evolutionary population dynamics with phenotypic fluctuation 
and derive an equation for the time evolution of population distribution in 
genotypic space. In Sec 13. H and Sec 13.21 we introduce a simple periodic fitness 
landscape as a multi-peaked fitness landscape. In Sec. 13.11 we investigate 
whether a population can keep its descendants around a peak of fitness, and 
can sustain high average fitness. In Sec. 13.21 we investigate how phenotype 
fluctuation affects evolutionary rate by using a transformation of the model 
equation into the Schrodinger equation and adopting WKB approximation. 
In Sec. 13.31 we investigate a quasi-periodic fitness landscape. In Sec. HI we 
include concluding remarks. 



2.1 Modulated fitness by phenotypic fiuctuation 

Let us begin with a discussion of how fitness landscapes are modified by phe- 
notypic fluctuation (also see [271 EO, EH])- Here, we analyze asexual popula- 
tion dynamics of a single-locus quantitative genetics model. Suppose that an 
individual organism has two types of a one-dimensional quantitative trait: a 
heritable trait of genotype g and a non-heritable trait of phenotype x. Geno- 
type g is inherited from its parent, whereas phenotype x is given through a 
genotype-phenotype mapping function p{x\g\ S) with fluctuation magnitude 
E. Fitness is assigned to each phenotype x and thus is expressed as F{x)^ 
given that natural selection acts on phenotypes x rather than genotypes g. 
By calculating the average with respect to x, the effective fitness assigned to 
the genotype g is given as the following equation: 



2 Model 




(1) 
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For simplicity, p{x\g; S) is assumed to be a Gaussian distribution with mean 
g and variance S; 

= ^==e . (2) 
From this simphfication, Eq. [1] leads to 

The above equation reveals that phenotypic fluctuations make a fltness land- 
scape flatter, and its peaks, lower. Examples of such modulation effects are 
shown in Fig. 1(a) and Fig. 1(b): By flattening the fltness landscape, fltness 
peak values are decreased, and thus it acts as a cost of phenotypic fluctuation. 

2.2 The Discrete Time Model and The Approximated 
Continuous Time Equation 

Once F{x) is given, the effective fltness assigned to individual organisms with 
genotype g is obtained in Eq. ([3]). We here consider the evolutionary pop- 
ulation dynamics of individual organisms. Suppose that N^{g; S) is the 
frequency distribution of the population at t-th generation, and each indi- 
vidual can produce descendants that survive in the next generation at a rate 
f{g)/{f)- Here, (■) represents an average over the entire population, (e.g., 
(/) = J^dgN*{g;T.)f{g;T.)/ J N^{g]T.)dg). Due to mutations, genotypes of 
t + 1-th generation g' deviate from the mother's genotype g as g' = g + C,, 
where ^ is a random number drawn from a Gaussian distribution with a mean 
zero and variance Dg. We deflne this Dg as the mutation rate. This effect of 
mutations on the distribution N{g) can be written using the mutation op- 
erator Mg^g, as Mg^g>[N{g')] = J N{g') ^ e-^3-g')V2D,^gf_ gy combining 

' ' Y 2n Dg 

both effects of natural selection and mutations, whole evolutionary dynamics 
is expressed as 

7V«teE) = AW|«^|. (4) 

With an assumption of sufficiently small Dg, this discrete time equation (jl]) is 
approximated by a continuous time equation governed by non-linear Fokker- 



5 



Planck-like equation: 

— dt — = — (J) — ^' + WW^^'' ^ ^ ' ^- ^ ^ 

In the above equation, N^{g; S) is replaced by N{t, g; S), which is a function 
of continuous time. From this approximation, we can obtain time evolution 
of a set of moment of N{t, g; S); for instance 

d{g)_ ^ {gf)-{g){f) ... 
dt if) 

d{g') _ {g'f)-{g'){f) 

dt if) 

3 Results 



Dr (7) 



3.1 Localized-extended transition in a periodic fitness 
landscape 

Here we investigate the model under a simple multi-peaked landscape, namely, 
a periodic one given as F{x) = 1 + cosaa;. Hence, the effective fitness land- 
scape is expressed as f{g) = 1 -|- exp(— q;^E/2) cosag' from Eq. ([3]). First, we 
investigate whether the population initially located at a peak of fitness can 
stay around the peak over generations to maintain high fitness. To examine 
this, we assume that at the initial state t = all individuals in a population 
are located at g = 0, and that the genotype distribution can be approximated 
by a Gaussian distribution with a mean value G{t) (= (g)) and variance Vg{t) 
(= ((yf^) — (g)"^) at any generation. This simplification allows us to calculate 
the time evolution of N{t,g] S) from Eq. ([6]) and Eq. ([7]), resulting in 



ag 



dt l+e-«^(S+V9)/2^Qg 



dVg_ _ g^Vf cos(«G)e-°'(^+'^a)/^ 
dt - i+e-"'(^+^9)/2 cos aG ^' 

The average population fitness is thus obtained as 



(8) 



(/) = (1 + cos(a^)e^— ) = 1 + cos(aG)e ^ — . (g) 



= [1 + E,?=i ij.{-§^rM„{g)]S{g-g') ^ {l + ^^)S{g-g') is used. See 
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At the initial state t = 0, we assume G(0) = and V^(0) = 0. Thus, Eq. ([S]) 
leads to dG/dt = (i.e., G(t) = 0) and 



,2^2g-a2(E+V'3)/2 



(10) 



dt 




The above equation has a fixed point solution when 



«2(^;)2e-"'(s+K,*)/2 



(11) 



^ 1 + g-a2(s+y;)/2 ' 
where V^* gives the minimum of dVg / dt as given by 

2 (2 + iy(2e-2-"'s/2)^ 



) 



(12) 



Here, W{-) is the Lambert W-function defined by the inverse function of 
f{W) = We^ . A finite Vg at equilibrium indicates that the population is 
localized around one of the peaks of fitness. When the inequality (fTT]) is 
not satisfied, Vg in Eq. ( JTOj) always increases in time, and thus, the variance 
of N(t,g; E) goes to infinity, leading to a broad distribution of N(t,g; S) in 
genotype space. This qualitative change in the variance Vg can be referred 
to as the localized-extended transition. 

Simpler representation of the condition Eq. ffTTl) is obtained by using 
the approximation 1 + exp[— " ^ 1^ which is reasonable around a 

transition point because variance Vg is sufficiently large near the transition 
point. It yields the condition for a finite variance at equilibrium, namely, the 
condition for the localized state as follows: 



To examine the validity of these arguments, we compare Eq. (ITT!) and Eq. (1131) 
with an individual-based simulation (i.e., direct numerical simulation of Eq. (|3j)). 
Figure 2 shows that both the estimates agree rather well with the results from 
the individual-based simulation. 

It should be noted that the Gaussian approximation fails when the pop- 
ulation is separated into two subgroups and localized at different peaks of 
fitness. Because of this failure, the variance of population estimated in the 
individual-based simulation does not coincide with the approximated value 




^ -a'^T./2-2 
1 



(13) 
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Vg'' calculated as a solution of Eq. ( TTUi) . To remedy this failure, we take 
account of only one peak with the largest subpopulation, as shown in [B], in 
which the procedure to determine the boundary between the localized and 
extended phases from the individual-based simulation is also described. 

The average population fitness is theoretically estimated as follows: by 
solving dVg/dt = in Eq. ffTU]) by using a quasi- Newton method, we compute 
Vg'^, and then the average population fitness is obtained by Eq. The 
estimated average fitness is shown in Fig. 3. We also estimate the average 
population fitness by using an alternative approximation method, a harmonic 
potential approximation, which is given below (see Eq. (!20l) and also Eq. ( I37I) 
in[D|. These results are shown in Fig. 3, and agree rather well with the results 
obtained from the individual-based simulation. 

Both approximated estimations and results from the individual-based 
simulation show that in the extended phase, the average population fit- 
ness takes a value that is as low as that without selection pressure, namely, 
(/) = 1. Such a drop in a high fitness value can be interpreted as "error 
catastrophe", i.e., the population fails to sustain high fitness due to high 
mutation rate or large phenotypic fluctuations. In contrast, the population 
in the localized phase is concentrated around the peak position of the fitness, 
and thus can sustain a high average fitness. 

3.2 Escape rate analysis in periodic fitness landscape 

Next we focus on the population dynamics of the localized phase. By the 
inherent stochastic nature of the evolutionary process, the population shows 
diffusion in genetic space, and maintains a high fitness value. Figure 4(a) 
shows examples of time series of mean genotype value [G] obtained in the 
individual-based simulations, whereas the time course of its mean square 
displacement [G^] is plotted in Fig. 4(b). In the figure, [G^] increases in 
proportion to t^, indicating that even in the localized phase the genotype 
in population performs random walks by jumping across a fitness valley. 
Through such random walks, a population can search for novel genotypes 
with a higher fitness value in genotype space without suffering from the error 
catastrophe. Thus, the diffusion coefficient of this random walk gives a mea- 
sure for the rate of evolutionary innovations, in other words, evolvability. The 
diffusion process in genotype space constitutes a sequence of jumps across 
fitness valleys. As Fig. 1(b) shows, phenotypic fluctuation tends to make fit- 
ness valleys shallower, and thus enhances the jump probability from a peak 
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per unit time F defined as /T^/^ N{t, g; T^)dg ~ e i.e., 



to a neighboring peak. Here, we analytically estimate the jump probability 

7v/a 
-tt/ck 

-tt/cm 



r = -llog N{t,g-T)dg, (14) 



with initial condition N{t = 0,g = 0;S) = 1. Note that the jump probabihty 
from a peak of fitness to a neighboring peak per unit time is one half of 
the escape rate in which the jumps to both neighboring peaks have to be 
considered. 

With a separation ansatz N{t,g;T,) = e^^(f){g;T,), Eq. 1^ becomes the 
following eigenvalue problem: 

mg) = ^^jy^m + ^^^JigHig), (i5) 

where (j){g\Y\) is rewritten as (f){g) and D is defined as D = Dg/2. By inter- 
preting (/) as an external control parameter R = (/), we obtain 

A0(^?) = (fig) - R)m + D—f{g)cl){g), (16) 
where A is given as 

A = R\. (17) 

By adopting variable as i/(^?) = ^ -j^dg' B.ndi,{y) = e^^y^^''^ ^(j){g{y)), 
where = ^ = a/ dD , the equation (fT6l) is transformed into 

-V{y) + d—U = Aij, (18) 
where 

n.)^ -(/(.(.)) -i^.q^-qf^£). m 

If we use tsc = —iht and m^c = h'^/2d, this agrees with the familiar form of 
Schrodinger equation. 

The Schrodinger form Eq. f lTSl) with WKB approximation [HI |35] allows 
us to estimate the jump probability between a peak of fitness and another 
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peak as an escape rate from a double well potential, assuming that the prob- 
ability distribution is located only in the left well at t = 0, as is illustrated 
in Fig. 5. Details of the derivation are given in [Dl 

Now, we consider the same fitness landscape in Sec. I3.lt f{g) = 1 + 
Xcosag and x = exp(— where x decreases as phenotypic fluctuation 
E increases. Using quasi equilibrium approximation (see details in[D]), the 
average population fitness R for this landscape can be derived as follows: 



(20) 

This evaluation of average fitness R at the quasi equilibrium agrees rather 
well with the results from the individual-based simulation, as is shown in 
Fig. 3 (the thick lines represent the evaluation in Eq. (1201) ), indicating that 
the quasi equilibrium approximation is valid. 

Using WKB approximation and R estimated above, we can evaluate the 
jump probability per unit time F as 



Rn \ y/D Jb V l + xcosa^f 

where 



e(y) = -^-(i + x)«^x + ^«V+^YT^^^-^(i-cosa,) 

(22) 

Note that g = b is a. point where the integrand becomes zero i.e., b satisfies 

X(l -cosafe) + 6(6) = (23) 
The effective diffusion coefficient De is given by 

The factor 2 in the above equation appears since the escape rate from a peak 
of fitness to the left neighbor peak has to be incorporated, as well as that to 
the right neighbor peak. Equations ( !2T|) and Eq. (12^ describe the relation 
between phenotypic fluctuation S and evolvability D^. 
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To examine the validity of the estimation of the escape rate 2r in Eq fl2ip . 
we compare Eq. fl2ip with that obtained numerically by the individual-based 
simulation of Eq. (jlj) (details of the estimation of the escape rate from the 
individual-based simulation are explained in[C]). As shown in Figs. 6(a) and 
6(b), the two agree rather well. Large deviation is discernible only for small 
S in Fig. 6(a), where steepness of potential V{x) invalidates the condition 
for WKB approximation. Relatively small deviation is also discernible for 
large E and large a in Fig. 6(a) and 6(b), since the double well potential 
approximation becomes less valid for lower effective fitness f{g) (i.e., shal- 
lower potential V{y) in Eq. (fT9|) ). These results shown in Figs. 6(a) and 6(b) 
indicate that phenotypic fluctuation always accelerates the evolutionary rate, 
while large phenotype fluctuation decreases the average population fitness, 
as is shown in Fig. 3. Therefore, the evolvability and fitness values are in a 
trade-off relationship. 



3.3 Mixed cosine potential 

To demonstrate that our analysis is not limited to the simple periodic po- 
tential, we consider the quasi-periodic fitness landscape; F{x) = Ai{l + 
cosaix) + ^2(1 + cosQ;2a;). From Eq. ([1]), the effective fitness landscape can 
be calculated as 

f{g, S) = Ai{l + cosaig) + ^2(1 + e 2" cos 025'). (25) 



An example of the effective fitness f{g) is illustrated in Fig. 7. 

Using Gaussian approximation of N{g; E), we obtain the time evolution 
of mean value G and variance Vg of Nig;!!) as 

~ ^Ai(l+e-°?(^+^9)/2cosaiG)+A2{l+e-"2(^+^f)/2cosa2G) /qcN 
~ 3 ^^(i+e-"^(^+^f)/2cosaiG)+A2(l+e-"2(^+^9)/2cosa2G) ^ 

Nullcline analysis of Eq. (1261) is given in Fig. 8(a) - (c), against the three 
values of S: for small S (Fig. 8(a)), a population can be localized at around 
each local maximum of fitness; for middle value of S (Fig. 8(b)), a population 
can be localized only at around a larger local maximum of fitness; for large 
S (Fig. 8(c)), all fixed points disappear, and Vg diverges, which indicates the 
extended phase. These results suggest the existence of two transitions: the 
localized-localized transition and the localized-extended transition. 
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(i) Localized-localized transition 

For the simplicity of discussion, we assume a\ > a2 in the following argument. 
Here, we define the localized-localized transition as the transition where a 
population cannot be localized at around the lowest local maximum of V{x). 
Thus, we investigate whether a population can be localized around Gim, 
which satisfies 

cos{aiGim) = 1, (27) 
cos{a2Gim) = -1. (28) 
In other words, we determine the conditions on which Vg has finite value for 
t — oo for a population at G = Gim from the following equation: 

dt ^9 ^j(i+g-<^f(S+V5)/2)_^^^(^_g-a|(E+V9)/2) ^9- y ^) 

Solving the above equation with an initial condition Vg{0) =0 by numerical 
integration, we determine the transition point where Vg diverges. In Fig. 9, 
the blue line shows the localized-localized transition point obtained by the 
numerical integration. 

Intuitive estimation is provided by assuming that the critical value of Vg 

is Vg ~ (207) ' ^^^'^'^ corresponds to the square of the width of a smaller 
peak. From this assumption, the transition point can be evaluated by 

' - ^2a,> Ai(l + e^"^(^+(^)^)/^) + ^(l-e-"^(^+(^)^)/^)' ^ ^ 
The above estimate agrees well with the results of the numerical integration 
of Eq. (129]) (see Fig. 9). 

(ii) Localized-extended transition 

To investigate the localized-extended transition, we consider the behaviors 
of the population distribution located at the global maximum G = 0. We as- 
sume that, near the transition point, Vg is large enough to consider e~°i'^^"'"^s^/^ 
as negligible, but small enough to interpret e~°2(^+^9)/2 ^ finite value. 
From this assumption, Eq. (!26l) can be rewritten as 

Here, the same procedures to derive Eq. ( fTTj) and Eq. ( fT3ll are applied. In 
Fig. 9, we show the estimated localized-extended transition point. 
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4 Discussion 



In the present study, the evolutionary population dynamics of a simple quan- 
titative genetics model under a multi-peaked fitness landscape is investigated 
by focusing on the role of phenotypic fluctuation on the evolutionary pro- 
cess. The main results of this study are summarized as the following points 



(i) The relationship between average fitness and phenotypic fluctuation. 
Analytical expression is obtained either by Gaussian approximation of popu- 
lation distribution or by a harmonic potential approximation. By the Gaus- 
sian approximation, average fitness is given by Eq. ([9]) with a solution of ODE 
Eq. ( fTOl) that requires numerical computation, while by the harmonic poten- 
tial approximation, the average fitness is analytically estimated as Eq. (1201) 
(or Eq. (JD])). As is shown in Fig. 3, both approximations agree rather well 
with the results from the individual-based simulation of evolutionary popu- 
lation dynamics described by Eq. and indicate that average population 
fitness decreases as the magnitude of phenotypic fluctuation S increases. 

(ii) Error catastrophe due to large phenotypic fluctuation. By either a 
high mutation rate or large phenotypic fluctuation, a population fails to keep 
its descendants around a peak of fitness, and thus is unable to sustain a high- 
fitness phenotype. Such a decrease in fitness is known as "error catastrophe." 
Although the error catastrophe usually refers to the failure in maintaining 
the population due to a high mutation rate, our result demonstrates that 
the error catastrophe can also occur by large phenotypic fluctuation. In our 
model, phenotypic fluctuation always enhances the risk of error catastrophe, 
whereas an interesting exception to this has been reported in which pheno- 
typic fluctuation can suppress the error catastrophe under a step-like fitness 
landscape in a high-dimensional genotype space |36j . 

(iii) Dependence of evolutionary rate on the magnitude of phenotypic 
fluctuation. We provide an analytical expression of the rate at which a pop- 
ulation jumps from one peak to another peak in the multi-peaked fitness 
landscape, which is expressed as the evolutionary rate. The obtained expres- 
sion Eq. fl2Tl) shows that phenotypic plasticity generally enhances evolution- 
ary rate. Such an acceleration has been noted previously [22], [271 1211 [25| [26]: 
however, our results provide the first analytical expression including both 
the effects of natural selection and mutation. Our claim here is that the 
acceleration of evolution by phenotypic fluctuation is different from Fisher's 
fundamental theorem [37j, which claims that genotypic variance is propor- 
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tional to evolutionary rate. In Fisher's theorem, the phenotypic plasticity or 
fluctuations are not considered, and the effect of mutation is also ignored, 
whereas in our results, genetic variance, phenotypic effect, and the effect of 
mutation are incorporated. A possible relationship between the phenotypic 
fluctuation and genotypic variance was also discussed in [5]. 

(iv) Phenotypic fluctuation exhibits a trade-off relationships with the ac- 
celeration of evolution and degradation of fitness. This remark is derived 
by combining the results (i), (ii), and (iii). This trade-off relationship is il- 
lustrated by both Fig. 3 and Fig. 6(a) (also see Fig. 6(b)), in which larger 
phenotypic fluctuation results in a smaller fitness value and a larger evolu- 
tionary rate. In addition, the transition to the extended phase (the error 
catastrophe) for larger phenotypic fluctuation can be interpreted as an ex- 
treme situation of this trade-off. It is commonly known that an increase in 
the mutation rate leads to the same trade-off relationship between the evolu- 
tionary rate and degradation of fitness, and that the error catastrophe occurs 
as an extreme case of the trade-off. Our results indicate that phenotypic fluc- 
tuation introduces an additional mutation rate, although the origin of the 
phenotypic fluctuation is essentially different from that of mutation. Note 
that this "additional mutation rate" effect of phenotypic fluctuation works 
only in case of multi-peaked fitness landscape, but not in case of a unimodal 
fitness landscape [2nlfT8l[2T]. 

Throughout the present study, we interpret random phenotype fluctua- 
tion as a kind of phenotypic plasticity (i.e., the phenotype distribution of 
a single genotype is independent of environment and fitness (see Eq. ([2]))). 
Several preceding studies addressing the Baldwin effect also adopted ran- 
dom fluctuation as phenotypic plasticity [271 1201 [HI [25] . Our results here 
indicate that even random phenotypic fluctuation accelerates the evolution- 
ary rate under a multi-peaked fitness landscape. Another interpretation 
of the phenotypic plasticity is that the phenotype distribution of a given 
genotype depends on the fitness landscape, referred to as responsive plas- 
ticity [ini [231 [221 [251 [26]. The difference between phenotypic fluctuation 
and responsive plasticity is only in the modulation of a fitness landscape 
by phenotypic changes, and therefore, the techniques that we adopted in 
the present study (e.g., the transformation of Fokker-Planck-like equation 
into the Schrodinger equation, WKB approximation) are applicable after a 
modulated fitness is obtained; these techniques are generally useful for inves- 
tigations of evolutionary dynamics [38l [39l [40l [36]. Although our analytical 
results are obtained for periodic and quasi-periodic fitness landscapes, they 
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would also be applicable to a general rugged landscape. 

In this paper we study an infinite population model in which population 
is always genetically polymorphic and there is no fixation event of a single 
genotype. It will be important for future research to extend our study to the 
case with a finite population, in which both the genetic drift and phenotypic 
fiuctuations have to be taken into account seriously. 

So far, mounting experimental evidences support that fitness landscapes 
in real organisms are multi-peaked H2l H3] . In addition, it has been 
also shown that stochasticity and random fiuctuations in phenotypes are 
common in living systems [Sj HI O [6] . Under the presence of such fiuctuation 
and multi-peaked fitness landscape, our analytical study has depicted the 
relationship between phenotypic fiuctuation and evolutionary rate, as well 
as the relationship between the Baldwin effect and the error catastrophe. 
Therefore, our findings have broad implications for the study of evolutionary 
dynamics. 
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Appendix 



A Detail of the individual-based simulation 

We use an individual-based evolutionary population dynamics model based 
on an algorithm similar to a genetic algorithm (GA); this simulation is a 
direct computation of Eq. 03]). In this model, each individual in a population 
has a genotype g and an effective fitness /((?), which is calculated in Eq. ([1]). 
We use = l + cos(ax) in Sec. 13.11 - Sec. l3.2l and = y4i(l + cosaix) + 
A2{1 + cosa2a;) in Sec. 13.31 

In our simulation, total population size is fixed to A^. At the first gen- 
eration, each individual has the same genotype g = Q. At the t-th gen- 
eration, A^ individuals for the next generation are selected from the cur- 
rent generation; each individual is sampled with probability f{g)/Nf, where 
/ = f{9i)/^- Mutations are applied to a selected population in a way 
in which g of each selected individual is mutated as g ^ g + ^. .^isa ran- 
dom number drawn from a Gaussian distribution /^'^^ / ^2i:Dg, where 
Dg is mutation rate. After a selection and mutation procedure, we obtain 
the t + 1-th generation. Iterating the procedure, we simulate the population 
dynamics. 



B Determination of the phase boundary in 



To restore the failure of the Gaussian approximation, we take account of only 
one interval with an integer n, Sn = {g\ {2nn — 7i)/a < g < (27m + n)/a} 
that has the largest subpopulation. At each time step in the simulation, 
variance of the population Vgim is evaluated by using the subpopulation in 
Sn by the following equation. 



Fig. 2 




(32) 



16 



Here, ns indicate the number of individuals in S'„. We also evaluated V* by 
truncated Gaussian distribution as 

^y* _ i:il9'eM-9V2V;)/,/2W;dg 

where V* is calculated by Eq. f|T2|) . In Fig. 2, localized (extended) phase is 
determined by the condition that the time average of Kim is smaller (larger) 
than V*. 

Note that the above estimate of Vsim using only subpopulation in un- 
derestimates the variance of the population, because broader distributions 
of coexisting subpopulations in neighboring peaks are ignored. The agree- 
ment between the boundaries by theory and simulation in Fig. 2 could be 
improved by taking appropriately into account a correction for the estimate 
of the variance. 



C Details of the evaluation of the escape rate 
from the individual-based simulation 

To evaluate the escape rate of a population from a peak of fitness, we per- 
form the individual-based simulation of Eq.(jl]) with the following boundary 
condition: at the initial condition t = 0, all individuals in a population are 
located at (7 = 0; when an individual escapes from the region (— tt/q;, tt/q;), 
g of the individual is transformed into g + Mir/a, where M is a sufficiently 
large integer that the individual never returns to the region {—Ti/a^iT/a). 
Note that this transformation does not change the fitness of each individ- 
ual. Using this boundary condition, we compute time series of the average 
fraction of a population in the region (— 7r/a,7r/a) and then estimate the 
decreasing rate of the average fraction as the escape rate. Here, the average 
fractions are calculated over 200 trials for a = 1.0,1.2 and 3000 trials for 
a = 0.8. Then the escape rates are calculated as 2r in Eq. flT^ . 



D WKB approximation 

Here, we assume that the population is located only around a peak of fitness. 
This assumption allows us to suppose that the probability distribution of 
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the Schrodinger equation is only located in the left well at t = 0, as is 
illustrated in Fig. 5. Within a time-scale much shorter than that of a jump, 
we can assume that a population is approximately in an equilibrium state (a 
quasi equilibrium) in a harmonic potential, which is derived from the Taylor 
expansion around a minimum oi V{y), y = y* 

V{y) ^ V{y*) + \v"{y*){y-y*f = V{y*) + \ + {y-yl\ 

(34) 

where S(?/) = f{y)— R. The second term in Eq. is interpreted as 
^rriscoo'^y'^, and therefore, 

where u is the angular frequency of the harmonic potential. The ground 
energy of the harmonic potential is given by 

(36) 

At the quasi equilibrium, the eigenvalue Eq should be zero, yielding 

-/V).q^-^|(-H..,.M) 

where /(y*) = f{g{y*)) = f{g*), E"iy*) = ^0|,=,* and U"{y*) = 
2 ag2 \g=9*- 

The escape rate of the double well potential can be calculated by the 
difference AA between the smallest eigenvalue Aq and the second smallest 
eigenvalue Ai [381 [32]. To derive this difference, the most successful technique 
is WKB approximation, in which sufficiently small and smooth potential 
V{ii) are assumed. By WKB approximation, the escape rate T of double well 
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potential can be estimated as [2S1 EH] 



AA hoj [ 2v2m 



exp 



2R R-K ' \ h 



^V{y) - Eody j = —exp (^--^ | ^^(y) - Eo 



(38) 

where |/ = 6 is a point in which the integrand becomes zero and y = a is a 
point in which the integrand becomes a maximum (i.e., the central peak of 
potential in Fig. 5 ). It should be noted that, from Eq. (ITTl) . the time scale 
of Eq. ([S]) is R times slower than Eq. fll6p . The integrand in Eq. f l38l) can be 

estimated as 



VV{y) -Eo = J- (siy) 



U"{y) {U'{y)y 

2 M 



Ad 



f{g*) - f{y) + ©(?/), (39) 



where 

<s>iy) = -1 



Dfig*)d^fi9*) , D\d^fig*)^^^ D .df{g)^,D.d'f{g*) d'f{g) 



2 ' 16' ^ 8/' 9^ ' 4' Qg2 

(40) 



Thus, equation Eq. (138|) becomes 



r = 1^ exp A v//(^7(r)) - /(^7(i/)) + e(2/)rf2/j 



v/2D/(r)i/"(r)l^^p ( _ J_ r n9:)^JM±^M9)ldg I 



^ ^ (41) 

Note that, in the above equation, we use huj = ^J D\f"{g*)\{2f{g*) + £)|/"(5(*)|/4) 

v/2D/((7*)|/"(^7*)| from Eq. ((351). 
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Figure Captions 



Fig. 1 Fitness landscapes, (a) A unimodal fitness landscape with (dot- 
ted line) and without (line) phenotypic fluctuation. Phenotypic fluctuation 
modulates an original fitness landscape to that with a lower peak and a wider 
range, (b) A multi-modal fitness landscape with (line) and without (dotted 
line) phenotypic fiuctuation. 

Fig. 2 A phase diagram for the localized-extended transition. 

The black line indicates the phase boundary estimated in Eq. flTTi) and the 
dotted line indicates the boundary estimated in Eq. ( fT3l) . The boundary 
between the blue and red regions is estimated from the individual-based sim- 
ulation. Details of the estimation are described in [Bl 

Fig. 3 Phenotypic fluctuation decreases average population flt- 
ness. The fitness are plotted against the phenotype fiuctuation S for Dg = 
0.2 and a = 0.8, 1, 1.2. Circles, squares and diamonds indicate results from 
the individual-based simulation. The dashed line indicates the approxi- 
mated value obtained from Eq. ([9]), where Vg'^ is calculated numerically from 
Eq. ( fTOj) . Thick lines indicate the approximated value obtained from a har- 
monic potential approximation in Eq. (l20l) . The deviation between dashed 
lines and thick lines will vanish for smaller Dg. Eq. (1131) indicates that the 
error catastrophe occurs as S > 4.76 for a = 1.0, as S > 2.80 for a = 1.2 
and as S > 8.84 for a = 0.8 . 



Fig. 4 Time series of the distribution of population in genotype 
space, (a) A time series of G (average g over population at each time) ob- 
tained by the individual-based simulation, (b) A log-log plot of the time 
series of [G^], where [] indicates the sample average. In both cases, the aver- 
age [■] is computed over 100 samples. [G"^] for each E can be fitted to a line 
with a slope equal to one. This indicates that G performs random walks. In 
both figures, = 500 and Dg = 0.3 are used. 
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Fig. 5 A double well potential. The probability is localized at the 
left well in the initial condition. A jump probability among peaks of fitness 
can be interpreted as the escape rate from the left well of the double well 
potential. 

Fig. 6 Escape rate from a peak of fitness, (a) The escape rate 2T for 
periodic potential F{g) = 1 + cos ag with respect to phenotype fluctuation 
E. Circles, squares and diamonds indicate results from the individual-based 
simulations. Dashed lines indicate the values estimated by using the WKB 
approximation in Eq. (|2T1) . (b) The escape rate 2T plotted against the phe- 
notype fluctuation E for the periodic potential F{g) = 2 + cos ag. Circles, 
squares and diamonds indicate results from the individual-based simulations. 
Dashed lines indicate values estimated by the WKB approximation. In both 
cases, N = 500 and Dg = 0.2 are used. 

Fig. 7 A mixed cosine fitness landscape. The landscpe is given by 
Eq. (!25|) . with (dotted line) and without (line) phenotypic fluctuation. Pa- 
rameters Ai = 2.5, A2 = 1.0, tti = 2.1, Q!2 = 1 are used. 

Fig. 8 Nullcline analysis for Eq. ([26]). Ai = 2.5, A2 = 1.0, ai = 
2.1,^2 = 1 and Dg = 0, 1 are used. The horizontal axis represents G, and 
the vertical axis, Vg. Black dots indicate stable fixed points, (a) for small E 
(E = 0), stable fixed points appear around each local maximum of fitness, 
(b) For E = 0.5, some of the stable fixed points vanish, and thus, there are 
only fixed points around G = ^ x n, n = (0, ±1, ±2, ...). This indicates that 
a population is localized only around a peak associated with 02. (c) For large 
E (E = 4), all fixed points vanish. In this case, the population goes to the 
extended phase. 

Fig. 9 A phase diagram of the localized-extended and localized- 
localized transition for a mixed cosine fitness landscape. Parameters 
are set as Ai = 2.5, A2 = 1.0, ai = 2.1, 02 = 1, and Dg = 0.1. The 
boundary between the localized phases 1 and 2 is estimated by numerical 
integration of Eq. ( l29l) (blue line) and by Eq. (1301) (red line). The boundary 
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between the localized phase 2 and the extended phase is estimated hy D = 

16 A2 -a\Y.l2-2 
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